TIFR/TH/03-07, lhep-lat/0303013 



Pressure and non-linear susceptibilities in QCD at finite chemical 

potentials 

Rajiv V. GavaQ and Sourendu GupteQ 
Department of Theoretical Physics, Tata Institute of Fundamental Research, 
Homi Bhabha Road, Mumbai 400005, India. 

Abstract 

When the free energy density of QCD is expanded in a Taylor series in the chemical potential, 
\i, the coefficients are the non-linear quark number susceptibilities. We show that these depend on 
the prescription for putting chemical potential on the lattice, making all extrapolations in chemical 
potential prescription dependent at finite lattice spacing. To put bounds on the prescription depen- 
dence, we investigate the magnitude of the nondinear susceptibilities over a range of temperature, 
T, in QCD with two degenerate flavours of light dynamical quarks at lattice spacing 1/4T. The 
prescription dependence is removed in quenched QCD through a continuum extrapolation, and the 
dependence of the pressure, P, on \i is obtained. 
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One of the most important objects in the study of hot and dense hadronic matter is 
the phase diagram, particularly, the location of the critical end point, characterised by the 
temperature T% and the chemical potential Much effort has been expended recently 
on estimating these quantities at finite lattice spacing, a, using, implicitly |l| or explicitly 
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a Taylor series expansion of the free energy density. This needs the non-linear 
susceptibilities which define the response to an applied \i beyond quadratic order. An 
equally important question for phenomenology arises from the fact that present day heavy- 
ion collision experiments access the part of the QCD phase diagram with \x ~10-80 MeV, 
i.e., baryon chemical potential jiB —30-250 MeV Q], far from \le- It is then pertinent to 
ask how relevant the \i = lattice QCD computations of quantities such as the pressure, P, 
are to these experiments. 

In this paper we present the first investigation of these non-linear susceptibilities. We 
uncover essential lattice artifacts, but manage to quantify and remove them in the process 
of taking the continuum limit. We explicitly construct a Taylor series expansion for P at 
H > 0, put limits on the region of linear response, i.e., of reliable extrapolations, and show 
that the /z = lattice computations are clearly relevant to experiments. An interesting 
sidelight is that there is strong evidence of short thermalisation times in the dense matter 
formed in these heavy- ion collisions which may be related to large values of transport 
coefficients 3]. Most computations of such dynamical quantities are based on linear response 
theory. The success of the linear approximation in static quantities at fairly large driving 
also gives us confidence in using linear response theory for dynamics. Another interesting 
point is that the radius of convergence of a Taylor series expansion started near T c |q] 
must give information on the location of the critical end-point, (T#, /i^), through the Taylor 
coefficients, i.e., the non-linear susceptibilities. Since these Taylor coefficients turn out to 
be prescription dependent and subject to strong finite lattice spacing effects, it seems that 
present day estimates of the end point will have to be sharpened strongly before they can 
be used as a guide to phenomenology. 

The partition function of QCD at finite temperature T and chemical potentials /i/ for 
each flavour / can be written as 

Z = e~ F/T = f VUe- s{T) Y[DetM(m f ,T,ij f ). (1) 

/ 

F is the free energy, S is the gluon part of the action, M is the Dirac operator, each 
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FIG. 1: All topologies which contribute to derivatives up to fourth order, and the notation for the 
corresponding operators. 



determinant is for one quark flavour and the temperature T enters through the shape of 
the lattice and boun dary conditions [9j. We shall work with a lattice discretisation and 
use staggered quarks [10(. In this work we shall only consider two degenerate flavours of 
quarks — m u = rrid = m [llj], with chemical potentials fi u and fid- The number densities, 
rif, and the (linear) quark number susceptibilities, Xfgi are the nrs t an d second derivatives 
of — F/V with respect to /// and \i g jl^. Since P = —F/V for a homogeneous system, the 
non-linear susceptibilities of order n > 3 are also the remaining Taylor coefficients of an 
expansion of P — 

1 d n F T d n \ogZ 
Xf 9 - = -77^7^7. = 777^7^ > ( 2 ) 



V dfijdfig ■ ■ ■ V dfMfdfjLg- • ■ 
where we construct the expansion around fif = 0. 

We now write systematic rules for the construction of the non-linear susceptibilities. The 
derivatives of log Z needed in eq. (J2J) can be related to the derivatives of Z with respect to 
the chemical potentials fjj, fi g , etc., (which we denote by Zf g ...) by the usual formulas for 
taking connected parts 13]. The only extra point to remember is that all the odd derivatives 
vanish by CP symmetry. To write the subsequent formulae compactly, we define operators 
0< by 



Z(0i), 



and 







rt+l 



dfif' 



(3) 



where angular brackets denote averages over the ensemble defined by eq. ((TJ) at /// = 0. 
Diagrammatic rules Q] for the 0j and the derivatives of Z, are- 

1. Put down n vertices (each corresponding to a derivative of M with respect to /i/) and 
label each with its flavour index. 
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2. Join the vertices by lines (each representing a quark) into sets of closed loops such 
that each loop contains only vertices of a single flavour. 0j is denoted by a single loop 
joining % vertices. 

3. For degenerate flavours and fif = 0, the operators are labeled only by the topology, 
which is specified completely by the number of vertices per loop and the number of 
such loops. Therefore erase the flavour index after step 2. We denote each resulting 



operator by the notation { 



0j0j • • •, where i + j + 



n. 



4. For each n-th order derivative of Z, add all the operator topologies for fixed n with 
flavour-dependent multiplicity equal to the number of ways in which each topology 
arises given the flavour indices. 



The number densities n u - 
(linear) susceptibilities X3 



series of papers 
derivatives 



rid = (T/V)(0i) vanish at n = 0. We have considered the 
- (T/V){0 2 ) and Xud = (T/V)(0n) extensively in a recent 



211 ] . The new quantities that we now consider are the two third order 



Z UUU = Z(0 3 + 30 12 + 1U ) and Z uud = Z(0 12 + 1U ), (4) 



the three fourth order derivatives 



Zuuuu = Z(04 + 40 13 + 30 22 + 60 11 2 + 0mi), 

Zuuud = ^(013 + 30H2 + 0llll), 
Zuudd = Z '(022 +20112 + 01111), 



(5) 



and the five corresponding susceptibilities. The third order susceptibilities turn out to 
vanish. The fourth order susceptibilities are 



Xuuuu 
Xuuud 
Xuudd 



z 

^uuud 

z 

^uudd 



z 



-3 



J ud 



(6) 



The operators contributing to eqs. (@J EJ) are shown in Figure H Note the interesting fact 
that beyond the second order, the number of distinct operator topologies is greater than the 
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tibilities Q; 



number of susceptibilities |14j |: however by making Nf sufficiently large, all topologies up to 
any given order can be given a physical meaning. 

A perturbative expansion in the continuum proceeds through an order-by-order enumer- 
ation of interaction terms. In the continuum the diagrams in Figure Q are the leading order 
(ideal quark gas) part of the perturbative expansion of the susceptibilities, where each ver- 
tex corresponds to the insertion of a 70 (since the chemical potential enters the Lagrangian 
as jofif)- Higher order Feynman diagrams correspond to dressing these loops by gluon 
attachments in all possible ways. 

In the lattice theory the diagrams in Figure d stand for operator definitions which need 
further specification. They are not Feynman diagrams, but mnemonics for the process of 
taking derivatives of Z. Since, the coupling of Fermions to the chemical potential is non- 
linear Q|, hence all derivatives of M exist and are non-zero in general. Using the identity 
DetM = exp(Tr InM) it is easy to get the usual expression &\ = Tr M _1 M', where M' is 
the first derivative of M with respect to a chemical potential. Next, using the chain rule 

d -^- = -M~ l M'M-\ (7) 
djif 

which comes from the identity MM -1 = 1, we recover the relation 2 = 
1 1 ( M L M'M ■ M' + M M"). where M" is the second derivative of M with respect 
to the chemical potential. Higher operators can be derived by repeated application of the 
chain rule with eq. ((ZJ), and involve higher derivatives of M, which we write as (a 
systematic method for doing this is given in the appendix). In particular, 



03 = Tr 

04 = Tr 



2(M~ 1 M') 3 - 3M l M' + M _1 M (3) 

- 6(M" 1 M') 4 - 3(M" 1 M") 2 + 12M- 1 M"(M~ 1 M') 2 
- \M ] M' :V: M 1 M' + M- l M^ 



(8) 



This completes the lattice definitions of the operators. 

Before we proceed to evaluate them and extract the non-linear susceptibilities, we note 
an ambiguity that arises on the lattice due to the fact that there is no unique way of putting 
chemical potential on the lattice. One can associate a factor f(a/i) for the propagation of a 
quark forward in time by one lattice spacing and a factor g(a/i) for the propagation of an 
antiquark. There are exactly four physical conditions that these two functions must satisfy 
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0. 



[151 ] . In the absence of chemical potential the usual lattice theory must be recovered, hence 
/(0) = g(0) = 1. CP symmetry gives /(— a/i) = <?(a/i). Finiteness of the energy density 
is guaranteed if f(a^i)g(afi) = 1. Finally, the correct continuum limit requires /'(0) = 1. 
These constraints imply the further relations, f"(0) = 1 and /^(0) = (— l) n g^ n \0), where 
the superscript n on / and g denotes the n-th derivative. All this guarantees that rif and 
Xf g are prescription independent. 

The four conditions above also give relations between the remaining f n \ such as = 
4j(3) _ 3^ ]3 U ^ q|q no t fix their numerical values. Since /i appears linearly in the continuum 
Lagrangian, these higher derivatives are all lattice artifacts. Any extra conditions imposed 
to fix them cannot be physical, and must remain at the level of prescription. The usual 
prescription, /(a/i) = exp(a/i) [16j . which we call the HK prescription, gives /^(0) = 1, 
but the alternative BG prescription /(a/i) = (1 + a/i)/-V(l — a 2 /i 2 ) 17| gives /^(0) = 3 and 
/ (4) (0) = 9. 

The difference between the two prescriptions can be rather significant. At any fixed cutoff, 
one may try to roughly map two prescriptions on to each other by changing /i while holding 
Z fixed by keeping /(a/i) unchanged. This gives the relation that for constant physics we 
must have 

a[i BG = tanh(a/iffif), (9) 
where this mapping is for quark chemical potentials. On N t = 4 lattices, the critical end- 



point for 2+1 flavour QCD has been determined to be at Te = 160 ± 3.5 MeV and /t^ 



725 ± 3 MeV 



HK 



i] 



The matching formula of eq. then shows that /if G ~ 692 MeV, and 
hence the resultant uncertainty in from this source alone is about 11 times larger than 
the statistical errors. We next show that this ambiguity vanishes in the continuum limit in 
all prescriptions. We also show later (Table EJ) that uncertainties of almost 20% are also 
expected from other finite lattice spacing effects even within one prescription, and lattice 
spacings of 1/12Te may be required to find fiE stable within statistical error bars. 

This freedom of choosing a prescription has specific consequences for the third and higher 
derivatives of M, and through them for the non-linear susceptibilities, and hence for F, P 
and all quantities at finite \x and a. At fif = 0, the derivatives of M are related by 

M (n) = f(n) a n-2 M , ^ odd ) M (n) = ^ eyen y (10) 

As a result, 3 = 0f K + Af^a 2 1 and 4 = 0f K + 4Af^a 2 2 , where the superscript 



HK on an operator denotes its value obtained in the HK prescription and A/® = /( 3 ) — 1. 
Clearly, the prescription dependence, manifested as a non-vanishing Af^ at this order, 
disappears in the continuum limit, a — > 0. Since (0!) = at \i = 0, the prescription 
dependence of (0 3 ) is invisible. We find that Xuuud = Xuuud + ^f {3) {Xud/T 2 )/N 2 . Since 
Xud vanishes within errors, as we show later, Xuuud turns out to be effectively prescription 
independent. From the relation for 04 we find, on varying N t at fixed T, 

y = \ HK + Af (3) (— ](— \ (11) 

Xuuuu A.uuuu ' J y rp2 J \j\f2 J ' \ 1A -J 

Finally, XuuM involves neither nor and hence is prescription independent. The 

prescription dependence of other susceptibilities can be systematically worked out, and it 
can be shown exactly as above that they become physical only in the continuum. Mixed 
derivatives of T and \x also have similar behaviour. If the dependence on a of each sus- 
ceptibility were known in any scheme, then one could write down an improved prescription 
by removing finite a effects systematically. In other schemes every quantity is potentially 
prescription dependent at finite lattice spacing. 

After this analysis of lattice artifacts in the Taylor coefficients, we return to the Taylor 
expansion itself. Along the line fi u — ^d — the Taylor series expansion of P can be written 
in the form 



AP fXuu\ (\i 



\H*/TJ V/< 



(12) 



2 / ,, ,"r \ 2 / , , r 

where AP = P(/i) — P(/i = 0), we have neglected Xuudd in anticipation of our numerical 
results ( Tables HI and ITT]) . and 



/i* _ 12x^/P 2 , , 

T ~ \ \y I ' 1 ' 

^ \ \Auuuu\ 

For an ideal gas in the continuum, Xuu/T 2 = 1 and Xuuuu = 6/7r 2 , giving /x*/T = \[2tx ~ 4.43. 
Some remarks are in order — 

1. The series within square brackets in eq. (|12|) is prescription dependent at any non-zero 
lattice spacing, and hence physical values of AP can be most reliably extracted by 
extrapolating each term in the series to the continuum. 

2. For those values of \i at which the second or higher order terms in the brackets in eq. 
(|12|) are important, computations of AP/T 4 on lattices with finite N t are necessarily 
prescription dependent. Since F = —PV, the same is evidently true for all other 



T/T c 


my/T 


W 6 Xud/T 2 


10 ^Cuuud 


10 ^Cuudd 


l£ K /T 


1.0 


0.2 


6 (30) 


4 (17) 


7(1) 


3.20 (3) 




0.1 


8 (42) 


7 (33) 


9(2) 


3.31 (5) 




0.03 


11 (84) 


20 (172) 


11 (2) 


3.38 (4) 


1.5 


0.2 


-0.3 (423) 


-0.7 (116) 


0.107 (3) 


3.73 (1) 




0.1 


0.6 (431) 


-0.6 (128) 


0.105 (3) 


3.84 (1) 




0.03 


-0.07 (433) 


-0.5 (166) 


0.106 (3) 


3.86 (2) 


2.0 


0.2 


2 (36) 


0.5 (85) 


0.097 (3) 


3.83 (1) 




0.1 


2 (36) 


0.5 (89) 


0.098 (3) 


3.87 (1) 




0.03 


1 (35) 


0.6 (82) 


0.096 (3) 


3.78 (2) 


3.0 


0.2 


0.6 (19) 


0.1 (5) 


0.032 (2) 


3.87 (1) 




0.1 


0.6 (20) 


0.1 (5) 


0.033 (2) 


3.88 (2) 




0.03 


0.6 (20) 


0.1 (5) 


0.033 (2) 


3.88 (2) 



TABLE I: Results in two flavour QCD with sea quark m/T c = 0.1. For T = T c the results are 
based on 2017 configurations, for 1.5T C on 370, for 2T C on 126 and for 3T C on 60. At T c and 3T C 
100 noise vectors were used. Xuuuu can be extracted from /x* and Xuu using eq. (fT3|) . 

physical quantities, including the energy density. From eqs. (fTD and IT^j) . it is clear 
that the prescription dependence of the quadratic term is (fi*/T) 2 . For N t — 4 
this can be as large as 33% (see Table EJ). 

3. If the series in eq. (II 2|) is well behaved, i.e., sixth and higher order susceptibilities are 
not much larger than Xuuuu, then this expansion must be well approximated by the 
leading term for u <C /i* in every prescription, and hence be effectively independent 
of prescription |l8| . Other finite lattice spacing effects may still exist. 

4. The series expansion must fail to converge in the vicinity of a phase transition; therefore 
estimates of (T E ,fi E ) on finite lattices must be prescription dependent, as we have 
already estimated. Computation of the continuum limit of several terms in the double 
series for F(T, /i) may allow us to use series extrapolation methods, such as Pade 
approximants or estimates of radius of convergence [3], to identify (T E ,[i E ) in the 
continuum limit. 
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W 6 Xud/T 2 
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IT 


1.5 


4 


2 (28) 


-0.7 (56) 


1.48 


(2) 


3.81 


(2) 




8 


0.2 (15) 


0.2 (13) 


0.70 


(1) 


4.36 


(4) 




10 


-0.4 (77) 


0.04 (64) 


0.61 


(2) 


4.47 


(4) 




12 


-0.5 (5) 


0.00 (30) 


0.56 


(1) 


4.55 


(4) 




14 


0.9 (58) 


0.00 (24) 


0.53 


(!) 


4.56 


(4) 
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0.45 


(!) 


4.67 


(4) 


2.0 


6 


0.2 (67) 


0.2 (10) 


1.01 


(1) 


4.11 


(1) 
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-0.3 (115) 


0.1 (13) 


0.74 


(1) 


4.32 


(5) 
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-0.3 (76) 


0.007 (49) 


6.37 


(3) 


4.45 


(3) 




12 


0.0 (57) 


0.00 (34) 


0.58 


(!) 


4.56 


(3) 




14 


-0.2 (43) 


0.00 (17) 


0.56 


(!) 


4.59 


(4) 
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0.49 


(3) 


4.76 


(4) 


3.0 


4 


2 (25) 


0.8 (44) 


1.54 


(1) 


3.85 


(1) 
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2 (4) 


0.1 (4) 


0.79 


(2) 


4.25 


(5) 




10 


-0.6 (14) 


-0.04 (11) 


0.66 


(!) 


4.40 


(3) 




12 


-0.1 (17) 


-0.02 (7) 


0.61 


(!) 


4.48 


(4) 




14 


-0.2 (8) 


0.00 (3) 


0.58 


(!) 


4.51 


(2) 




oo 






0.496 (1) 


4.62 


(1) 



TABLE II: Results in quenched QCD with m v /T c = 0.1. Quadratic extrapolations to the con- 
tinuum limit, Nt = oo, from the last three points, are shown, fi* and x,uuuu related by eq. 



We turn now to our numerical simulations. For dynamical sea quark mass m/T c = 0.1 we 
studied the higher order susceptibilities at T = T c on a 4x 10 3 lattice, 1.5T C and 2T C on 4x 12 3 
lattices, and 3Tf on a 4 x 14 3 lattice. All the simulations were performed using the hybrid 
R-algorithm 20j with molecular dynamics trajectories integrated over one unit of MD time 
using a leap-frog algorithm with time step of 0.01 units. At T c autocorrelations of the Wilson 
line and the quark condensate were found to be between 150 and 250 trajectories. With 
over 2000 saved configurations separated by 10 trajectories each, this gave the equivalent 
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of about 100 independent configurations. For T > T c the autocorrelations were all less 
than 10 trajectories, and hence all the saved configurations can be considered statistically 
independent. 

Quark number susceptibilities were evaluated in the HK prescription on stored configu- 
rations using valence quark masses my/T c = 0.2, 0.1 and 0.03. The smallest valence quark 
mass is chosen such that the ratio of the (T = 0) rho and pion masses reaches its physical 
value 0.2 at the lattice spacing a = 1/4T C . All quark-line disconnected diagrams of the 
kind needed for these measurements are evaluated using a straightforward extension of the 



stochastic method given for Xud in 21] using 10 to 100 noise vectors per configuration 
Our results for the non-linear susceptibilities which do not vanish by symmetry are shown 
in Table E It is clear that of these only Xuudd and Xuuuu, are non-zero with statistical sig- 
nificance. Comparing them to computations with sea quark mass m/T c = 0.2 and various 
volumes, we concluded that they are free of sea quark mass and finite volume effects. Also 
note the stability in physical quantities as m v /T c decreases from 0.1 to 0.03. 

With present day computer resources the continuum limit is hard to take in QCD with dy- 
namical quarks. To investigate this limit we have evaluated the same quantities in quenched 
QCD for T > 1.5T C where the difference in the order of transitions is immaterial |23|. The 



run parameters are exactly as in 



2l| . Our results are shown in Table HTl These results show 



that there is over 20% movement in /i* when going from N t = 4 to the continuum within 
a fixed prescription. Since /i* is an estimate of the radius of convergence of the Taylor 
expansion at the fourth order, it implies that the estimate of the end-point, /jle, may shift 
upward by about 20% due to finite size effects even inside the HK scheme. Xuudd remains 
significantly non-zero on all the lattices, and there is some evidence that it becomes either 



zero or marginally negative in the continuum 



24( | . We shall present more detailed studies in 



the future. Finally, the results for N t = 4 are very similar in the quenched and dynamical 
theories, leading us to believe that the continuum limits will also be close. 

^ obtained in inched QCD, U8 i„ g vaMes of tan Q and ,,/T obtained 
here, are shown in Figure El At RHIC it is seen that fi/T c = 0.06 <C 0.15, which implies 
that AP/T 4 is negligible. In terms of dimensionless variables, the results in quenched and 

23- For fi/T c ~ 0.45, 



dynamical QCD are not expected to differ by more than 5-10% 
relevant to SPS energies, the effects of fi > are more significant, but can still be reliably 
extracted using only the leading term of eq. (112(1 . In this whole range of [ijT c the results 
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FIG. 2: AP/T 4 as a function of T/T c for the values of fi/T c shown. Continuum results correct 
to 0(/^ 4 ) (full lines) and 0(/U 2 ) (dotted lines) are shown. Nt = 4 results, in the HK prescription, 
correct to 0(^ 4 ) and multiplied by 0.47 to compensate for finite a effects in Xuu are shown with 
dashed lines. 

of [25!, including a correction for finite lattice spacing artifacts in the evaluation of Xuu at 
Nt = 4, are the same as our continuum results, and both are dominated by the leading 
term of eq. (|12p. Our computations show that for fi > 2T C , higher order terms become 
significant for the continuum limit. As a result, at these chemical potentials, reweighting on 
Nt = 4 lattices, even after correcting for finite a effects in Xuu, are quite different from the 
continuum values. 

In conclusion, we have studied non-linear susceptibilities and shown that they are pre- 
scription dependent at finite lattice spacing. We have found the continuum limit of these 
quantities in quenched QCD, and thereby removed these artifacts. This allows us to compute 
the finite chemical potential corrections to the pressure relevant to RHIC and SPS experi- 
ments. For a = 1/4T the numerical results for QCD with and without dynamical quarks are 
similar, and we find the continuum limit of some of these quantities in the quenched theory. 
It would be interesting to compare them with perturbation theory. We have argued that 
the critical end point (Te,He) evaluated at N t = 4 is uncertain by more than 10 times the 
statistical errors. As a result, a continuum extrapolation is required to obtain the physical 
value of the end point. This may be possible with the computation of several non-linear 
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susceptibilities. 

We would like to thank J.-P. Blaizot for discussions. 



APPENDIX A: LATTICE OPERATORS 

In this appendix we work in lattice units, i.e., we choose the lattice spacing to be unity. 
We introduce the compact notation 



Tr 



m -i m (pi) M^M^ 



(nrpi©n2-pj©--), (Al) 



and further write (1 ■ p) as (p). Since the trace allows only cyclic permutations, therefore 

(affi&ffic) = (cffiaffifr) ^ (frffiaffic), (A2) 

i.e., the 'addition' (represented by ©) is not commutative. 'Multiplication' (denoted by the 
dot) is distributive over addition, subject to restrictions due to non-commutativity, i.e., 

(n-p@m-p) = ((n + m)-p), 
in • p © m ■ p' @l • p) = ((n + I) ■ p © m ■ p'), (A3) 

but no simplification is possible for (n • p © m • p' ffi I • p © • • •). Traces can be added, i.e., 

a(n ■ p) + b(n ■ p) = (a + b)(n ■ p). (A4) 

The point of all this is to simplify the taking of derivatives. These are easy to write — 

(n . p y = - n (\ © n ■ p) + n((n - 1) • p © (p + 1)). (A5) 

The operation of taking derivatives is linear over the 'addition' ©, since this is just the rule 
for taking derivatives of products. 
We have the first examples 

1 = (1), 2 = -(2-1) + (2). (A6) 

Then, the remaining known ones are obtained simply by applying the rules again. Since 
(2-1)' = -2(3- 1)+ 2(1 ©2) and (2)' = -(1 ffi 2) + (3), we first obtain the relation in eq. ©, 

3 = 2(3-1) -3(1 ffi 2) + (3). (A7) 
12 



At the fourth order we need the derivatives 



(3 - 1)' = -3(4- 1) + 3(2-1 ©2), 

(1 © 2)' = -2(2 • 1 © 2) + (2 • 2) + (1 © 3), 

(3)' = -(10 3) + (4). (A8) 

As a consequence of the general rule in eq. (jA5|) . the coefficients sum up to zero. This 
is a consequence of the rule for derivatives in eq. (|A5j) . Also note that each operator, 
(• • • © rij • pi © • • •), which contributes to n must satisfy the constraint S^-iPi = n. The 
expressions in eq. (jA8|) give the result of eq. (jHJ) , 

4 = -6(4 ■ 1) + 12(2 ■ 1 © 2) - 3(2 ■ 2) - 4(1 © 3) + (4). (A9) 

For each n for n > 2, the sum of the coefficients is zero, as can be proved by induction 
from eq. (|A5|) . 

Using these rules higher order derivatives, needed for the higher order susceptibilities, can 
be easily written down. Since these manipulations are simple rules for rewriting expressions, 
not only are they easy to automate inside standard algebra packages, but can even be readily 
implemented as macros in text editors such as sed or emacs. 
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